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A grid-free variant of the Direct Simulation Monte Carlo (DSMC) method is proposed, 
named the Isotropic DSMC (I-DSMC) method, that is suitable for simulating dense fluid 
flows at molecular scales. The I-DSMC algorithm eliminates all grid artifacts from the 
traditional DSMC algorithm; it is Galilean invariant and microscopically isotropic. The 
stochastic collision rules in I-DSMC are modified to yield a non-ideal structure factor that 
gives consistent compressibility, as first proposed in [Phys. Rev. Lett. 101:075902 (2008)]. 
The resulting Stochastic Hard Sphere Dynamics (SHSD) fluid is empirically shown to be 
thermodynamically identical to a deterministic Hamiltonian system of penetrable spheres 
interacting with a linear core pair potential, well-described by the hypernetted chain (HNC) 
approximation. We apply a stochastic Enskog kinetic theory to the SHSD fluid to obtain 
estimates for the transport coefficients that are in excellent agreement with particle sim- 
ulations over a wide range of densities and collision rates. The fluctuating hydrodynamic 
behavior of the SHSD fluid is verified by comparing its dynamic structure factor against 
theory based on the Landau-Lifshitz Navier-Stokes equations. We also study the Brownian 
motion of a nano-particle suspended in an SHSD fluid and find a long-time power-law tail 
in its velocity autocorrelation function consistent with hydrodynamic theory and molecular 
dynamics calculations. 

With the increased interest in nano- and micro-fluidics, it has become necessary to develop 
tools for hydrodynamic calculations at the atomistic scale [TJ [2]. There are several issues present 
in microscopic flows that are difficult to account for in models relying on the continuum Navier- 
Stokes equations. Firstly, it is complicated to deal with boundaries and interfaces in a way that 
consistently accounts for the bidirectional coupling between the flow and (moving) complex surfaces 
or suspended particles. Furthermore, it is not trivial to include thermal fluctuations in Navier- 
Stokes solvers [31 E| |5], and in fact, most of the time the fluctuations are not included even though 
they can be very important at instabilities [6] or in driving the dynamics of suspended objects 
[3 E]. Finally, since the grid cell sizes needed to resolve complex microscopic flows are small, a 
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large computational effort is needed even for continuum solvers. An alternative is to use particle- 
based methods, which are explicit and unconditionally stable and rather simple to implement. The 
fluid particles are directly coupled to the microgeometry, for example, they directly interact with 
the beads of a polymer chain. Fluctuations occur naturally and the algorithm may be designed to 
give the correct spatio-temporal correlations. 

Several particle methods have been described in the literature. The most accurate but also most 
expensive is molecular dynamics (MD) [9], and several coarse-grained models have been developed, 
such as dissipative particle dynamics (DPD) [10] and multi-particle collision dynamics (MPCD) 
|11| , I12|. each of which has its own advantages and disadvantages [13 . Our method, first proposed 
in Ref. [13], is based on the Direct Simulation Monte Carlo (DSMC) algorithm of Bird [15]. The 
key idea behind DSMC is to replace deterministic interactions between the particles with stochastic 
momentum exchange (collisions) between nearby particles. While DSMC is usually viewed as a 
kinetic Monte Carlo algorithm for solving the Boltzmann equation for a low-density gas, it can also 
be viewed as an alternative to the expensive MD in cases where an approximate (coarse-grained) 
treatment of the molecular transport is appropriate. The stochastic treatment of collisions makes 
the algorithm much simpler and faster than MD, while preserving the essential ingredients of 
fluctuating hydrodynamics: local momentum conservation, linear momentum exchange on length 
scales comparable to the particle size, and a similar fluctuation spectrum. 

Being composed of point particles, the DSMC fluid has no internal structure, has an ideal gas 
equation of state (EOS), and is thus very compressible. As a consequence, the density fluctuations 
in DSMC are significantly larger than those in realistic liquids. Furthermore, the speed of sound 
is small (comparable to the average speed of the particles) and thus subsonic (Mach number less 
than one) flows are limited to relatively small Reynolds numbers 1 . Efforts have been undertaken to 
develop coarse-grained models that have greater computational efficiency than brute-force MD and 
that have a non-ideal EOS, such as the Lattice-Boltzmann (LB) method [TB], DPD [T7], MPCD 
\18\ I19j. The Consistent Boltzmann Algorithm (CBA) |20[ 121], as well as algorithms based on the 
Enskog equation |22|. 123] . have demonstrated that DSMC fluids can have dense-fluid compressibility, 
however, they did not achieve thermodynamic consistency between the equation of state and the 
fluid structure. 

For a low-density gas the Reynolds number is Re « M/K, where M = Vfi om /c is the Mach number, and the 
Knudsen number K — X/L is the ratio between the mean free path A and the typical obstacle length L. This 
shows that subsonic flows can only achieve high Re flows for small Knudsen numbers, i.e., large numbers of DSMC 
particles. 
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In this paper we describe a generalization of the traditional DSMC algorithm suitable for dense 
fluid flows. By a dense fluid we mean a fluid where the mean free path is small compared to the 
typical inter-atomic distance. As a first step, we introduce a grid-free Isotropic DSMC (I-DSMC) 
method that eliminates all grid artifacts from traditional DSMC, notably the lack of Galilean 
invariance and non-isotropy. The I-DSMC fluid is still an ideal fluid just like the traditional 
DSMC fluid, that is, it has an the equation of state of an ideal gas and does not have an internal 
structure as do liquids. Secondly, by biasing the collision kernel in I-DSMC to only allow stochastic 
collisions between approaching particles, we obtain the Stochastic Hard-Sphere Dynamics (SHSD) 
algorithm that is thermodynamically consistent (i.e., the direct calculation of compressibility from 
density fluctuations agrees with the density derivative of pressure) . The SHSD algorithm is related 
to previous algorithms for solving the Enskog kinetic equation |22| [23] . and can be viewed as a 
more-efficient variable-diameter stochastic modification of the traditional hard-sphere molecular 
dynamics |24] . 

In the SHSD algorithm randomly chosen pairs of approaching particles that lie less than a 
given diameter of each other undergo collisions as if they were hard spheres of diameter equal to 
their actual separation. The SHSD fluid is shown to be non-ideal, with structure and equation 
of state equivalent to that of a deterministic (Hamiltonian) fluid where penetrable spheres effec- 
tively interact with a repulsive linear core pairwise potential. We theoretically demonstrate this 
correspondence at low densities. Remarkably, we numerically find that this effective interaction 
potential, similar to the quadratic core potential used in many DPD variants, is valid at all den- 
sities. Therefore, the SHSD fluid, as DPD, is intrinsically thermodynamically-consistent since it 
satisfies the virial theorem. 

The equivalence of the structure of the SHSD fluid with the linear core fluid enables us to use 
the Hypernetted Chain (HNC) approximation, as recommended in Ref. [25], to obtain theoretical 
estimates for the pair correlation and static structure factor that are in excellent agreement with 
numerical results. These further enable us to use the Enskog- like kinetic theory developed in Ref. 
|26j to obtain accurate theoretical estimates of the transport properties of the SHSD fluid that 
are also shown to be in excellent agreement with numerics even at relatively high densities. At 
lower densities the HNC approximation is not necessary and explicit expressions for the transport 
coefficients can be obtained similarly to what has been done using Green-Kubo approach for other 
DSMC variants [20, and MPCD [HKIHIET]. 

We numerically demonstrate that the hydrodynamics of the SHSD fluid is consistent with the 
equations of fluctuating hydrodynamics when the appropriate equation of state is taken into ac- 
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count. Specifically, we compare the measured dynamic structure factors with that obtained from 
the linearized fluctuating Navier-Stokes equations. We also calculate the velocity autocorrelation 
function (VACF) for a large hard spherical bead suspended in an SHSD fluid, demonstrating the 
existence of long-time tails as predicted by hydrodynamics and found in MD simulations. The tail 
is found to be in quantitative agreement with theory at lower densities, but a discrepancy is found 
at higher densities, possibly due to the strong structuring of the dense SHSD fluid. 

We begin by introducing a grid-free variant of the DSMC algorithm in Section [I] This Isotropic 
DSMC algorithm simulates a stochastic particle system where particles closer than a particle di- 
ameter collide with a certain rate. By biasing the collision kernels to favor head-on collisions of 
particles, as in the hard-sphere fluid, we obtain a non-ideal stochastic fluid in Section [TTJ We de- 
velop an Enskog-like kinetic theory for this Stochastic Hard Sphere Dynamics (SHSD) system in 
Section II A, which requires as input the pair correlation function. In Section II B| we discover that 
the SHSD fluid is thermodynamically consistent with a fluid of penetrable linear core spheres, and 
use this equivalence to compute the pair correlation function of the SHSD fluid using the HNC 



approximation. In Section III we show several numerical results, including a comparison with the- 
ory for the transport coefficients and for the dynamic structure factor, as well as a study of the 
hydrodynamic tails in the velocity autocorrelation of a bead suspended in an SHSD fluid. 



I. ISOTROPIC DSMC 



The traditional DSMC algorithm |15[ [28] starts with a time step where particles are propa- 
gated advectively, r i = + V{At, and sorted into a grid of cells. Then, for each cell c a certain 
number N co u ~ T SC N C {N C — l)At of stochastic collisions are executed between pairs of particles 
randomly chosen from the N c particles inside the cell, where the collision rate T sc is chosen based 
on kinetic theory. The conservative stochastic collisions exchange momentum and energy between 
two particles i and j that is not correlated with the actual positions of the particles. Typically the 
probability of collision is made proportional to the magnitude of the relative velocity v r = \vij\ by 
using a conventional rejection procedure. 

Traditional DSMC suffers from several grid artifacts, which become pronounced when the mean 
free path becomes comparable to the DSMC cell size. Firstly, the method is not Galilean invariant 
unless the grid of cells is shifted randomly before each collision step, as typically done in the MPCD 
algorithm [T2"] for the same reason. This shifting is trivial in a purely particle simulation with 
periodic boundary conditions, but it causes implementation difficulties when boundaries are present 
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and also in particle-continuum hybrids [29]. Furthermore, traditional DSMC, unlike MD, is not 
microscopically isotropic and does not conserve angular momentum, leading to an anisotropic 
collisional stress tensor. Instead of trying to work around these grid artifacts, as done for non-ideal 
MPCD in Refs. |18[ [19], we have chosen to modify the traditional DSMC algorithm to make the 
dynamics grid-free. 

To ensure isotropy, all particle pairs within a collision diameter D (i.e., overlapping particles if 
we consider the particles to be spheres of diameter D) are considered as potential collision partners 
even if they are in neighboring cells. In this way, the grid is only used as a tool to find neighboring 
particles efficiently, but does not otherwise affect the properties of the resulting stochastic fluid. 
Such a grid- free DSMC variant, which we will call the Isotropic Direct Simulation Monte Carlo 
(I-DSMC) method, is suitable for hydrodynamics of dense fluids, where the mean free path is 
comparable or even smaller than D, unlike the original DSMC which targets the dilute limit. It is 
important to point out, however, that the I-DSMC is not meant to be a replacement for traditional 
DSMC for rarified gas flows. In particular, the computational efficiency is reduced by a factor 
of 2 — 4 over traditional DSMC due to the need to search neighboring cells for collision partners 
in addition to the current cell. This added cost is not justified at low densities, where the grid 
artifacts of traditional DSMC are small. Furthermore, the I-DSMC method is not intended as a 
solver for the Boltzmann equation, which was the primary purpose of traditional DSMC |30l I3T] . 
Rather, in the limit of small time steps, the I-DSMC method simulates the following stochastic 
particle system: Particles move ballistically in-between collisions. While two particles i and j are 
less than a diameter D apart, nj < D, there is a probability rate xD -1 K c (vij,rij) for them to 
collide and change velocities without changing their positions, where K c is some function of the 
relative position and velocity of the pair, and the dimensionless cross-section factor x se ts the 
collisional frequency. Because the particles are penetrable, D and x may be interpreted as the 
range and strength, respectively, of the interaction potential. After the collision, the pair center- 
of-mass velocity does not change, ensuring momentum conservation, while the relative velocity 



\vij\\ so kinetic energy is 



is drawn from a probability density P^v^Vij^ij), such that 
conserved. 

Once the pre- and post-collision kernels K c and P c are specified, the properties of the resulting 
I-DSMC fluid are determined by the cross-section factor x an d the density (hard-sphere volume 
fraction) <fi = irND 3 /(6V), where N is the total number of particles in the simulation volume 
V. Compare this to the deterministic hard-sphere fluid, whose properties are determined by the 
volume fraction 6 alone. It is convenient to normalize the collision kernel K r so that for an ideal 
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gas with a Maxwell-Boltzmann velocity distribution the average collisional rate would be x times 
larger than that of a gas of hard spheres of diameter D at low densities, <ft <C 1. Two particular 
choices for the pre-collision kernel K c that we use in practice are: 

Traditional DSMC collisions (Traditional I-DSMC ideal fluid), for which the probability of colli- 
sion is made proportional to the magnitude of the relative velocity v re \ = ||«y||, K c = 3v re i/4. 
We use this kernel mainly for comparison with traditional DSMC. 

Maxwell collisions (Maxwell I-DSMC ideal fluid), for which K c = 3v re i/4 = SWksl^/TTrri, where 
v re i is the average relative velocity at equilibrium temperature To. Since K c is a constant, all 
pairs collide at the same rate, independent of their relative velocity. This kernel is not realistic 
and may lead to unphysical results in cases where there are large density and temperature 
gradients, however, it is computationally most efficient since there is no rejection based 
on relative velocity. We therefore prefer this kernel for problems where the temperature 
dependence of the transport properties is not important, and what we will typically mean 
when we say I-DSMC without further qualification. 

Other collision kernels may be used in I-DSMC, though we will not consider them here [28 . 
We typically chose the traditional DSMC post-collisional kernel P c in which the direction of the 
post-collisional relative velocity is randomized so as to mimic the average distribution of collision 
impact parameters in a low-density hard-sphere gas. Specifically, in three dimensions the relative 
velocity is rotated uniformly independent of [15]. If one wishes to microscopically conserve 
angular momentum in I-DSMC then the post-collisional kernel has to use the actual positions of 
the colliding particles. Specifically, the component of the relative velocity perpendicular to the line 
joining the colliding particles should remain unchanged, while the parallel component should be 
reversed. 

Note that a pairwise Anderson thermostat proposed within the context of MD/DPD by Lowe [32] 
adds I-DSMC-like collisions to ordinary MD. In addition to algorithmic differences with I-DSMC, 
in Lowe's method the post-collisional kernel is such that it preserves the normal component of the 
relative velocity (thus conserving angular momentum) , while the parallel component is thermalized 
by drawing from a Maxwell-Boltzmann distribution. We strive to preserve exact conservation of 
both momentum and energy in the collision kernels we use, without artificial energy transport via 
thermostating. 

With a finite time step, the I-DSMC method can be viewed as a time-driven kinetic Monte Carlo 
algorithm to solve the Master Equation for the stochastic particle system described above. Unlike 
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the singular kernel in the Boltzmann equation, this Master Equation has a mollified collision kernel 
with a finite compact support D |26l [33] . The traditional DSMC method also mollifies the collision 
kernel by considering particles within the same collision cell, of size L c , as possible collision partners. 
This DSMC cell size is much larger than a molecular diameter, D m , in fact, for low densities it is a 
fraction (typically a quarter) of the mean free path. The molecular properties enter in traditional 
DSMC only in the form of collisional cross-sections a ~ D^. In light of this, for rarified gas flows, 
the collision diameter D in I-DSMC should be considered the equivalent of the cell length L c , and 
not D m . Traditional DSMC is designed to reproduce a collision rate per particle per unit time 
equal to the Boltzmann rate, r#(L> m ) = CD^, where C is a constant. The I-DSMC method is 
designed to reproduce a collision rate 



we get r i-dsmc = T_b(-Dot)- Therefore, if I-DSMC is used to simulate the transport in a low- 
density gas of hard-sphere of diameter D m , the collision diameter D should be chosen to be some 
fraction of the mean free path A (say, D ~ A/4 3> D m ), and the cross-section factor set to 
Xb ~ (-Dm/A) 2 <C 1. At higher densities \B starts becoming comparable to unity and thus it is no 
longer possible to separate the kinetic and collisional time scales as assumed in traditional DSMC. 
Note that I-DSMC is designed for dense fluids so while it is possible to apply it in simulating 
rarefied gases it will not be as computationally efficient as traditional DSMC. 



In I-DSMC, stochastic collisions are processed at the begining of every time step of duration At, 
and then each particle i is streamed advectively with constant velocity Uj. During the collision step, 
we need to randomly and without bias choose pairs of overlapping particles for collision, given the 
current configuration of the system. This can be done, as in traditional DSMC, using a rejection 
Monte Carlo technique. Specifically, we need to choose a large number Nf c ot ^ = N pa i rs At of 
trial collision pairs, and then accept the fraction of them that are actually overlapping as collision 
candidates. Here N pa i rs is the number of possibly-overlapping pairs, for example, as a first guess 
one can include all pairs, N pa i rs = N(N — l)/2. The probability for choosing one of the overlapping 



r i-dsmc = x^b(d) = X CD 2 , 



and therefore by choosing 




A. Performing Stochastic Collisions 
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pairs as a collision candidate is simply Tj c 'At. If the probability of accepting a candidate pair 
ij for an actual collision is and At is sufficiently small, then the probability rate to actually 

collide particles i and j while they are overlapping approaches = p^^vf"^. The goal is to 
choose the trial collision frequency T^"^ and pf^ 1 such that = xD^ 1 K c (vij, r^). 

The efficiency of the algorithm is increased if the probability of accepting trial collisions is 
increased. In order to increase the acceptance probability, one should reduce N pa i rs to be closer to 
the number of actually overlapping pairs, ideally, one would build a list of all the overlapping pairs 
(making N pa i rs linear instead of quadratic in N). This is however expensive, and a reasonable 
compromise is to use collision cells similarly to what is done in classical DSMC and also MD 
algorithms. Namely, the spatial domain of the simulation is divided into cells of length L c > D, 
and for each cell a linked list C c of all the particles in that cell is maintained. All pairs of particles 
that reside in the same or neighboring cells are considered as potential collision partners, and here 
we include the cell itself in its list of neighboring cells, i.e., each cell has 3 d neighbors, where d is 
the spatial dimension. 

To avoid any spatial correlations (inhomogeneity and non-isotropy) , trial collision pairs should 
be chosen at random one by one. This would require first choosing a pair of neighboring cells with 
the correct probability, and then choosing a particle from each cell (rejecting self-collisions). This 
is rather expensive to do, especially at lower when few actual collisions occur at each time step, 
and we have therefore chosen to use a method that introduces a small bias each time step, but is 
unbiased over many time steps. Specifically, we visit the cells one by one and for each cell c we 

(c) (c) 

perform N^ c = T\ c N c N p At trial collisions between one of the iV c particles in that cell and one 
of the N p particles in the 3 neighboring cells, rejecting self-collisions. Here TJ C ; is a local trial 
collision rate and it may depend on the particular cell c under consideration. Note that each of 
the N C (N C — 1) trial pairs ij where both i and j are in cell c is counted twice, and similarly, any 
pair where % and j are in different cells c and d is included as a trial pair twice, once when each of 
the cells c and d is considered. Also note that it is important not to visit the cells in a fixed order 
during every time step. Unlike in traditional cells, where cells are independent of each other and 
can be visited in an arbitrary order (even in parallel), in I-DSMC it is necessary to ensure isotropy 
by visiting the cells in a random order, different every time step. 

For the Maxwell pre-collision kernel, once a pair of overlapping particles i and j is found 
a collision is performed without additional rejection, therefore, we set = x^~ 1 K c /2 = 

3x(2.D) _1 \J 'ksTo I 'irm = const; note that we have divided by two because of the double count- 
ing of each trial pair. For the traditional pre-collision kernel, and, as we shall see shortly, the 
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SHSD pre-collision kernel, additional rejection based on the relative velocity Vij is necessary. As 
in the traditional DSMC algorithm, we estimate an upper bound for the maximal value of the 
pre-collision kernel j{^ max ^ among the pairs under consideration and set T^s = X-D -1 .?^™^ /2. 
We then perform an actual collision for the trial pair ij with probability 

giving the correct collision probability for every overlapping pair of particles. For the traditional 
pre-collision kernel ^ max ) = 3t^ ax ^/4, where v^ ax ^ is as tight an estimate of the maximum 
relative speed as possible. In the traditional DSMC algorithm uji^ is a global bound obtained 
by simply keeping track of the maximum particle speed v max and taking v^" 3 ^ = 2v m ax PS] • In 
I-DSMC, we obtain a local estimate of i& for each cell c that is visited, thus increasing the 
acceptance rate and improving efficiency 

Algorithm [T] specifies the procedure for performing collisions in the I-DSMC method. The 
algorithm is to a large degree collision-kernel independent, and in particular, the same algorithm 
is used for ideal and non-ideal stochastic fluids. As already explained, the size of the cells should 
be chosen to be as close as possible but still larger than the particle diameter D. The time step 
should be chosen such that a typical particle travels a distance l&t ~ v^/St ~ D5t, where the 
typical thermal velocity vth = \/kBTo/m and 5t is a dimensionless time step, which should be 
kept reasonably smaller than one, for example, St is 0.25. It is also important to ensure that each 
particle does not, on average, undergo more than one collision per time step; we usually keep the 
number of collisions per particle per time step less than one half. Since a typical value of the 
pre-collision kernel is K c ~ vth, the number of collisions per particle per time-step can easily be 
seen to be on the order of 

Naps ~ X VJ ^ ■ yV p ■ At = X (f>St, 

where V p ~ D 3 is a particle volume. Therefore, unless \§ 3> 1, choosing a small dimensionless time 
step bt will ensure that the collisional frequency is not too large, N cps <C 1. With these conditions 
observed, we find little dependence of the fluid properties on the actual value of 5t. 

Algorithm 1: Processing of stochastic collisions between overlapping particles at a time-step in 

the I-DSMC method. 

1. Sample a random permutation of the cell numbering V . 

2. Visit the cells one by one in the random order given by V. For each cell c, do the following steps if 
the number of particles in that cell iV c > 0, otherwise move on to the next cell. 
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3. Build a list C\ of the N c particles in the cell and at the same time find the largest particle speed in 
that cell v™ ax . Also keep track of the second largest speed in that cell u™ *, which is an estimate of 
the largest possible speed of a collision partner for the particle with speed v™ ax . 

4. Build a list of the N p particles in the set of 3 d cells that neighbor c, including the cell c itself and 
respecting the proper boundary conditions. Also update v™ ax if any of the potential collision partners 
not in cell c have speeds greater than v™ ax . 

5. Determine the number of trial collisions between a particle in cell c and a neighboring particle by 
rounding to an integer |15j the expected value 

N tc = T tc N c N p At, 

where At is the time step. Here the local trial collision rate is 

_ X Kr x 

tc ~ 2D ' 

where K™ ax is an upper bound for the pre-collision kernel among all candidate pairs. For Maxwell 
collisions K™ ax — Z^fk^T^Jwm, and for traditional collisions K™ ax = 3u^ oa5 ' /4, where Vmax = 
(v™ ax + v™ ax ) is a local upper bound on the relative speed of a colliding pair. 

6. Perform trial collisions by randomly selecting N tc pairs of particles i € C± and j S £2- For each pair, 
do the following steps if i =/= j: 

(a) Calculate the distance Uj between the centroids of particles i and j, and go to the next pair if 
kj > D. 

(b) Calculate the collision kernel Kf^ = K c (vij, r^), and go to the next pair if Kfj = 0. 

(c) Sample a random uniform variate < r < 1 and go to the next pair if Kfj < rK™ ax (note that 
this step can be skipped in Maxwell I-DSMC since = K™ ax ). 

(d) Process a stochastic collision between the two particles by updating the particle velocities 
by sampling the post-collision kernel P c {v i j\Vij^rij). For ideal fluids we perform the usual 
stochastic DSMC collision by randomly rotating Vij to obtain t?-., independent of ry. 



II. STOCHASTIC HARD SPHERE DYNAMICS 

The traditional DSMC fluid has no internal structure so it has an ideal gas equation of state 
(EOS), p = PV/NksT = 1, and is thus very compressible. As for the classical hard-sphere fluid, 
the pressure of fluids with stochastic collisions consists of two parts, the usual kinetic contribution 
that gives the ideal-gas pressure Pk = 1, and a collisional contribution proportional to the virial 
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Pc ~ (( v ij ' r ij)' ~ ( v ij ' r ij)) c i where the average is over stochastic collisions and primes denote 
post-collisional values. The virial vanishes for collision kernels where velocity updates and positions 
are uncorrelated, as in traditional DSMC, leaving only the ideal-gas kinetic contribution. In order 
to introduce a non-trivial equation of state it is necessary to either give an additional displacement 
Arij to the particles that is parallel to Vij , or to bias the momentum exchange Ap^- = mAvij to be 
(statistically) aligned to r^. The former approach has already been investigated in the Consistent 
Boltzmann Algorithm (CBA) |20[ [21]. This algorithm was named "consistent" because both the 
transport coefficients and the equation of state are consistent with those of a hard-sphere fluid to 
lowest order in density, unlike traditional DSMC which only matches the transport coefficients. 
However, CBA is not thermodynamically consistent since it modifies the compressibility without 
affecting the density fluctuations (i.e., the structure of the fluid is still that of a perfect gas). 

Here we explore the option of biasing the stochastic momentum exchange based on the position 
of the colliding particles. What we are trying to emulate through this bias is an effective repulsion 
between overlapping particles. This repulsion will be maximized if we make Ap^ parallel to J*y, 
that is, if we use the hard-sphere collision rule P c (u ij -; Vij, r^) = 6(vij+2v n rij), where v n = —Vij-rij 
is the normal component of the relative velocity. Explicitly, we collide particles as if they are elastic 
hard spheres of diameter equal to the distance between them at the time of the collision, 

v\ =Vi + Vnfij 

Vj =Vj -VnVij. (1) 

Such collisions produce a positive virial only if the particles are approaching each other, i.e., if 
v n > 0, therefore, we reject collisions among particles that are moving apart, K c (vij,rij) ~ Q(v n ), 
where O denotes the Heaviside function. Note that the hard-sphere post-collision rule ([I]) strictly 
conserves angular momentum in addition to linear momentum and energy and can be used with 
other pre-collision kernels (e.g., Maxwell) if one wishes to conserve angular momentum. 

To avoid rejection of candidate collision pairs and thus make the algorithm most efficient, it 
would be best if the pre-collision kernel K c is independent of the relative velocity as for Maxwell 
collisions. However, without rejection based on the normal v n or relative v r speeds, fluctuations of 
the local temperature T c would not be consistently coupled to the local pressure. Namely, without 
rejection the local collisional frequency T sc would be independent of T c and thus the collisional 
contribution to the pressure p c ~ (AiJjj • ry) ~ r sc \/Ic would be p c ~ \fT c instead of p c ~ T c , 
as is required for a fluid with no internal energy [TT?] . Instead, as for hard spheres, we require 
that T sc ~ y/Tc, which is satisfied if the collision kernel is linear in the magnitude of the relative 
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velocity. For DSMC the collisional rules can be manipulated arbitrarily to obtain the desired 
transport coefficients, however, for non-ideal fluids thermodynamic requirements eliminate some 
of the freedom. This important observation has not been taken into account in other algorithms 
that randomize hard-sphere molecular dynamics |34j . but has been used in the non-ideal MPCD 
algorithm in order to obtain thermodynamic consistency [18, .19] • 

There are two obvious choices for a pre-collision kernel that are linear in the magnitude of the 
relative velocity. One is to use the relative speed, K c ~ v r , as in the traditional DSMC algorithm, 
and the other is to use the hard-sphere pre-collision kernel, K c ~ v n . We have chosen to make the 
collision probability linear in the normal speed v n , specifically, we take K c = 3v n Q(v n ) to define the 
Stochastic Hard-Sphere Dynamics (SHSD) fluid, similarly to what has previously been done in the 
Enskog DSMC algorithm [22j 123] and in non-ideal MPCD pS Ql|. These choices for the collision 
kernels make the SHSD fluid identical to the one proposed in Ref. [33] for the purposes of proving 
convergence of a microscopic model to the Navier-Stokes equations. Specifically, the singular 
Boltzmann hard-sphere collision kernel is mollified in Ref. [33] to obtain the SHSD collision kernel 
and then the low-density hydrodynamic limit is considered. 

The non-ideal SHSD fluid is simulated by the I-DSMC method, in the limit of sufficiently small 
time steps. However, it is important to observe that the SHSD fluid is defined independently of 
any temporal discretization used in computer simulations, just like a Hamiltonian fluid is defined 
through the equations of motion independently of Molecular Dynamics (MD). To summarize, in 
the SHSD algorithm we use the following collision kernels in Algorithm [T] 

K c =3v n Q(v n ) and = 3v£*> 

Pc(v' i:j ) =S(vij + 2v n rij) 
where v n = —Vij ■ rij. 

Note that considering particles in neighboring cells as collision partners is essential in SHSD in 
order to ensure isotropy of the collisional (non-ideal) component of the pressure tensor. It is also 
important to traverse the cells in random order when processing collisions, as well as to ensure 
a sufficiently small time step is used to faithfully simulate the SHSD fluid. Note that the SHSD 
algorithm strictly conserves both momentum and energy independent of the time step. 
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A. Enskog Kinetic Theory 

In this section we develop some kinetic equations for the SHSD fluid that are inspired by the 
Enskog theory of hard-sphere fluids. Remarkably, it turns out that these sorts of kinetic equations 
have already been studied in the literature for purely theoretical purposes. 



1. BBGKY Hierarchy 

The full Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy of Master equations de- 
scribing the SHSD fluid is derived in Ref. [33] . Specifically, the evolution of the s-particle distri- 
bution function f s (t; n, V\, . . . , r s , v s ) is governed by 



dt 



s fl r r s 

+ y^Vi-V r J s =2>xD 2 I dx I dvj dr™ x 2S ^v n 



[f s+ i(t;n,vi, ...,ri,Vi,.. .,r 8 ,v a ,ri + xrij,Vj) 
-f 8+1 (t;ri,vx, ...,r i} Vi,.. .,r 8 ,v a ,ri - xrij,Vj)] (2) 

which takes into account the contribution from collisions of one of the s particles, particle i, with 
another particle j that is at a distance = xD away, < x < 1. Here denotes the fraction 
of the unit sphere for which v n = —r^ ■ (v, — Vj) > 0, and v\ = V{ + v n Tij and v'j = Vj — v n rij. 
Just like the BBGKY hierarchy for Hamiltonian fluids, Eqs. ^ are exact, however, they form an 
infinite unclosed system in which the (s+ l)-particle distribution function appears in the equation 
for the s-particle distribution function. As usual, we need to make an anzatz to truncate and close 
the system, as we do next. 



2. Thermodynamic and Transport Properties 

The hydrodynamics of the SHSD fluid is well-described by a kinetic equation for the single- 
particle probability distribution f(t,r,v) = fi(t;r,v) obtained by making the common molecular 
chaos assumption about the two-particle distribution function, 

h(M ri,v 1 ,r 2 , v 2 ) = g2(ri,r 2 ; n)f(t, r 1 ,v 1 )f(t, r 2 , v 2 ), 

where g 2 (ri, r,-; n) is the non-equilibrium pair distribution function that is a functional of the local 
number density n(r). At global equilibrium n(r) = const and g 2 = g 2 {rij) depends only on the 
radial distance once the equilibrium density n and cross-section factor x are specified. Substituting 
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the above assumption for f'2 in the first equation of the BBGKY hierarchy we get a stochastic 
revised Enskog equation of the form studied in Ref. [26 , 



df{t,r,v) v r /(t,r,«) =3 X D 2 [ dx [ dw [ de x 2 v n 

Jo Jn 3 J si 



Of 

[g 2 (r, r + xe; n)f(t, r, v')f(t, r + xe, w') 
-52 (r, r - xe; n)f(t, r, v)f(t, r - xe, w)] (3) 

where v n = —e ■ (v — w) > 0, v' = v + ev n and w' = w — ev n . 

The standard second-order Chapman-Enskog expansion has been carried out for the "stochastic 
Enskog" equation of the same form as Eq. ([3]) in Ref. [26 , giving the equation of state (EOS) 
p = PV/NksT, and estimates of the diffusion coefficient £, the shear r] and bulk ijb viscosities, 
and thermal conductivity k of the SHSD fluid. The expressions in Ref. [26] ultimately express 
the transport coefficients in terms of various dimensionless integer moments of the pair correlation 
function g2(x = r/D), x^ = J Q x k g2(x)dx, specifically, 

p - 1 =120xa?3i (4) 



VB/Vo= §75—, (6) 

, 5 24</>xx 3 2 3 

= 7F^ v 1 + ? + i»7B> and ( 7 

48^/ttxx2 5 5 

25 36</>xx 3 2 3 

K/K0 =6V^ (1 + } + 2^' (8) 

where Co = D yJksT /m, i]o = D^yJmk^T and /to = kgD~ 2 \j 'k^I '/ 'm are natural units. These 

equations are very similar to the ones in the Enskog theory of the hard-sphere fluid except that 

various coefficients are replaced with moments of g2{x). In order to use these equations, however, 

we need to have a good approximation to the pair correlation function, i.e., to the structure of the 

SHSD fluid. It is important to point out that Eq. Q is exact as it can be derived directly from 

the definition of the collisional contribution to the pressure. 



B. Pair Correlation Function 



In this section we study the structure of the SHSD fluid, theoretically at low densities, and then 
numerically at higher densities. We find, surprisingly, that there is a thermodynamic correspon- 
dence between the stochastic SHSD fluid and a deterministic penetrable-sphere fluid. 
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1. Low Densities 

In order to understand properties of the SHSD fluid as a function of the density <f> and the cross- 
section factor Xi we first consider the equilibrium pair correlation function g 2 (r) at low densities, 
where correlations higher than pairwise can be ignored. We consider the cloud of point walkers 
ij representing the N(N — l)/2 pairs of particles, each at position r = r\ — r» and with velocity 
v = Vi — Vj. If one of these walkers is closer than D to the origin, r < D, and is approaching 
the origin, v n > 0, it reverses its radial speed as a stochastic process with a time-dependent rate 
r = \v n \ To, where Tq = 3x/D is the collision frequency. A given walker corresponding to pair ij 
also undergoes stochastic spatially-unbiased velocity changes with some rate due to the collisions 
of i with other particles. At low densities we can assume that these additional collisions merely 
thermalize the velocities to a Maxwell-Boltzmann distribution but not otherwise couple with the 
radial dependence of the one-particle distribution function f pa irs(v, r) of the N{N — l)/2 walkers. 
Inside the core r < D this distribution of pair walkers satisfies a kinetic equation 



I pairs " J pairs 



dt dr 



I"* J~pairs ^f *^ 

— ^O^nfpairsi (9) 

I" j pairs if V n < 



where the term —Tf pa i rs is a loss term for approaching pairs due to their collisions, while Tf pa i rs 
is a gain term for pairs that are moving part due to collisions of approaching pairs that then 
reverse their radial speed. At equilibrium, df pa i rs /dt = and v n cancels on both sides, con- 
sistent with choosing collision probability linear in \v n \, giving df pa irs/dr = 3x-D _1 fpairs- At 
equilibrium, the distribution of the point walkers in phase space ought to be of the separable form 
f P airs(v,r) = f P airs(v n ,r) ~ g 2 (r) exp(-mv%/ikT), giving dg 2 (r)/dr = 3xD~ 1 g 2 (r) for r < D and 
zero otherwise, with solution 



exp [3y(x — 1)1 for x < 1 

11 n ~ (10) 
1 for x > 1 



g 2 (x = r/D) -- 

Indeed, numerical experiments confirmed that at sufficiently low densities the equilibrium g 2 for 



the SHSD fluid has the exponential form (10) inside the collision core. From statistical mechanics 
we know that for a deterministic Hamiltonian particle system with a pairwise potential U (r) at low 
density g% = exp[— U(r)/kT]. Therefore, the low density result (10) is consistent with an effective 
linear core pair potential 

U eff (r)/kT = 3 X (1 - x)6(l - x). (11) 
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Note that this repulsive potential is similar to the quadratic core potential used in DPD and strictly 
vanishes outside of the overlap region, as expected. Also note that the cross-section factor x plays 
the role that U(0)/kT plays in the system of penetrable spheres interacting with a linear core 
pairwise potential. 

As pointed out earlier, Eq. Q is exact. At the same time, it is equivalent to the virial theorem 
for the linear core potential. Therefore, if the pair correlation functions of the SHSD fluid and the 
linear core fluid are truly identical, the pressure of the SHSD fluid is identically equal to that of the 
corresponding penetrable sphere system. As a consequence, thermodynamic consistency between 
the structure [52 (aO and S(k)] and equation-of-state \p(4>)] is guaranteed to be exact for the SHSD 
fluid. 



2. Equivalence to the Linear Core Penetrable Sphere System 



Remarkably, we find numerically that the effective potential (11 ) can predict exactly g2(x) at all 
densities. In fact, we have numerically observed that the SHSD fluid behaves thermodynamically 
identically to a system of penetrable spheres interacting with a linear core pairwise potential for 
all (f> and x- Figure [T] shows a comparison between the pair correlation function of the SHSD fluid 
on one hand, and a Monte Carlo calculation using the linear core pair potential on the other, at 
several densities. Also shown is a numerical solution to the hypernetted chain (HNC) integral 
equations for the linear core system, inspired by its success for the Gaussian core model [25]. The 
excellent agreement at all densities permits the use of the HNC result in practical applications, 
notably the calculation of the transport coefficients via the Enskog-like kinetic theory presented 



in Section II A 2 We also show the static structure factor S(k) in Fig. [T] and find very good 



agreement between numerical results and the HNC theory, as expected since S(k) can be expressed 
as the Fourier transform of h(r) = 52 (^) — 1- 

For collision frequencies x 1 t ne structure of the SHSD fluid is quite different from that of 
the hard-sphere fluid because the particles inter-penetrate and overlap significantly. Interestingly, 
in the limit of infinite collision frequency x ~* 00 t ne SHSD fluid reduces to the hard-sphere (HS) 
fluid for sufficiently low densities. In fact, if the density <p is smaller than the freezing point for the 
HS system, the structure of the SHSD fluid approaches, as x increases, that of the HS fluid. For 
higher densities, if x is sufficiently high, crystallization is observed in SHSD, either to the usual 
hard- sphere crystals if <f) is lower than the close-packing density, or if not, to an unusual partially 
ordered state with multiple occupancy per site, typical of weakly repulsive potentials [35]. Monte 
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Figure 1: (Left) Equilibrium pair correlation function of the SHSD fluid (solid symbols, N — 10 4 particles 
in a cubic periodic box), compared to Monte Carlo simulations (open symbols, N = 10 4 particles in a cubic 
periodic box) and numerical solution of the HNC equations (solid lines) for the linear core system, at various 



densities and x — 1- The low-density approximation corresponding to Eq. (10 1 is also shown. (Right) The 
corresponding static structure factors from SHSD simulations (solid symbols, average of ten snapshots of a 
system with N = 10 5 particles in a cubic periodic box) and HNC calculations (solid lines). The time step 
was kept sufficiently small in the SHSD simulations to ensure that the results are faithfully represent the 
SHSD fluid with time-step errors smaller than the statistical uncertainty. 



Carlo simulations of the linear core penetrable sphere system show identical freezing behavior 
with SHSD, confirming the surprising equivalence even at non-fluid densities. This points to a 
conjecture that the (unique) stationary solution to the BBGKY hierarchy Q is the equilibrium 
Gibbs distribution, 



/, 



E _ Itt=i M{vi) 
Zn 



r s+ i Jr N 



exp 



dr s+1 . . . dr N , 



where M is a Maxwellian. 



III. RESULTS 



In this Section we perform several numerical experiments with the SHSD algorithm. Firstly, we 
compare the theoretical predictions for the transport properties of the SHSD fluid based on the HNC 
theory for the linear core penetrable sphere system with results from particle simulations. We then 
compute dynamic structure factors and compare them to predictions of fluctuating hydrodynamics. 
Finally, we study the motion of a Brownian bead suspended in an SHSD fluid. 
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A. Transport Coefficients 



The equation of state of the SHSD fluid for a given x is P = p(<ft)NkBT/V , where p(<f>) is given 
in Eq. Q. According to statistical mechanics, the structure factor at the origin is equal to the 
isothermal compressibility, that is, 

S = S(u = 0, k = 0) = c^ 2 = (p + <j)dp/d<j))- 1 

where = cx^/ksT /m is the isothermal speed of sound. In the inset in the top part of Fig. [2J we 
directly demonstrate the thermodynamic consistency of SHSD by comparing the compressibility 
calculated through numerical differentiation of the pressure, to the structure factor at the origin. 
The pressure is easily measured in the particle simulations by keeping track of the total collisional 
momentum exchange during a long period, and its derivative was obtained by numerical differen- 
tiation. The structure factor is obtained through a temporal average of a Fast Fourier Transform 
approximation to the discrete Fourier Transform of the particle positions \\^2 i exp(—ik • n) || 2 . The 
value S(k = 0) is estimated by fitting a parabolic dependence for small k and extrapolating to 
k = 0. 







St=0.025 fluct. 


□ 


§t=0.1 fluct. 




51=0.25 fluct. 
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8t=0.S fluct. 




5t=Q.Q25 compress. 




5l=0.1 compress. 




5l=0.25 compress. 




5l=0.5 compress. 




Figure 2: (Left) Normalized equation of state (p — l)J(x<f>) for the SHSD fluid at several cross-section 
factors x (different symbols, N — 10 5 particles in a cubic box) compared to theoretical predictions based 
on the virial theorem Q with the HNC approximation to (x) (solid lines). The inset compares the 
compressibility (p + 4>dp / d<f))~ l (dashed lines) to the structure factor at the origin S(k — > 0) (symbols), 
measured using a direct Fourier transform of the particle positions for small k and extrapolating to k = 0. 
The dimensionless time step St — 0.025 is kept constant and small as the density is changed. (Right) 
Thermodynamic consistency between the compressibility (lines) and the large-scale density fluctuations 
S(k — > 0) (symbols) for different dimensionless time steps St, keeping x = 1 Axed. 



As pointed out earlier, the dimensionless time step St = D/^kBTo/m in the SHSD algorithm 
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should be kept reasonably small, St <C min [l, ((j>x)~\i i n order to faithfully simulate the SHSD 
fluid. As the time step becomes too large we expect to see deviations from the correspondence with 
the linear core system and thus a violation of thermodynamic consistency This is indeed observed 
in our numerical results, shown in Fig. [2j where we compare the structure factor at the origin as 
estimated through the equation of state with that obtained from a direct Fourier transform of the 
particle positions. We should point out that when discussing thermodynamic consistency one has 
to define what is meant by the the derivative dp/ dep. We choose to keep the collisional frequency 
prefactor x an d the dimensionless time step St constant as we change the density, that is, we study 
the thermodynamic consistency of a time- discrete SHSD fluid defined by the parameters 4>, x an d 
5t. The results in Fig. [2]show that there are significant deviations from thermodynamic consistency 
when the average number of collisions per particle per time step is larger than one. This happens 
at the highest densities for 5t > 0.1, but is not a problem at the lowest densities. Nevertheless, 
a visible inconsistency is observed even at the lower densities for St > 0.25, which comes because 
particles travel too far compared to their own size during a time step. 




Figure 3: Comparison between numerical results for SHSD at several collision frequencies (different symbols) 
with predictions based on the stochastic Enskog equation using the HNC approximation for gi (x) (solid 
lines). The low-density approximations are also indicated (dashed lines). (Left) The normalized shear 
viscosity r//?yo a t high and low densities (inset), as measured using an externally-forced Poiseuille flow. 
There are significant corrections (Knudsen regime) for large mean free paths (i.e., at low densities and low 
collision rates). (Right) The normalized diffusion coefficient C(x < / ) )/Co, a s measured from the mean square 
displacement of the particles. The time step was kept sufficiently small in the SHSD simulations to ensure 
that the results arc faithfully represent the SHSD fluid with time-step errors smaller than the statistical and 
measurement errors. 



Having established that the HNC closure provides an excellent approximation g\ ~ g2 
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for the pair correlation function of the SHSD fluid, we can obtain estimates for the transport 
coefficients by calculating the first four moments of g^ N ( x ) an d substituting them in the results 



of the Enskog kinetic theory presented in Section II A2[ In Figure [3] we compare the theoretical 
predictions for the diffusion coefficient £ and the viscosity r] to the ones directly calculated from 
SHSD particle simulations. We measure £ directly from the average mean square displacement 
of the particles. We estimate rj by calculating the mean flow rate in Poiseuille parabolic flow 
between two thermal hard walls due to an applied constant force on the particles 2 . Surprisingly, 
good agreement is found for the shear viscosity at all densities. Similar matching was observed for 
the thermal conductivity k. The corresponding results for the diffusion coefficient show significant 
(~ 25%) deviations for the self-diffusion coefficient at higher densities because of larger corrections 
due to higher-order correlations. 



B. Dynamic Structure Factors 



The hydrodynamics of the spontaneous thermal fluctuations in the SHSD fluid is expected 
to be described by the Landau-Lifshitz Navier-Stokes (LLNS) equations for the fluctuating field 
U = (Po + 3p,8v,To + ST) linearized around a reference equilibrium state Uq = (poi^o = 0,To) 
[3"S| 157]. For the SHSD fluid the linearized equation of state is 

p = PW) , r ~ iPo + c T — +po—)p c , 

v PO -tO 

and there is no internal energy contribution to the energy density, 

„ 3 Nk B T 
e «: - — — — = e + c v lodp + p c v dl, 

where po = p(4>o), Co = ksT/m, and c v = 3kB/2rn, giving an adiabatic speed of sound c s = c s cq, 
where c? a = + 2p 2 /3. Omitting the S's for notational simplicity, for one-dimensional flows the 
LLNS equations take the form 



dtp 




d t v 


d 


dx 


d t T _ 





p v 

4po 1 p + Po4 t o 1t 



PqcIc v 1 v 



+ 



d_ 
dx 



Po 1 Vov x 

-1 -1 m 

Po c v K oT x 



+ 



d_ 

dx 



(12) 



2 Similar results are obtained by calculating the viscous contributions to the kinetic and collisional stress tensor 
in non-equilibrium simulations of Couette shear flow. This kind of calculation additionally gives the split in the 
viscosity between kinetic and collisional contributions. 
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where l^W and are independent spatio-temporal white noise Gaussian fields. 

By solving these equations in the Fourier wavevector-frequency domain for U(k,u) and per- 
forming an ensemble average over the fluctuating stresses we can obtain the equilibrium (sta- 
tionary) spatio-temporal correlations (covariance) of the fluctuating fields. We express these cor- 
relations in terms of the 3x3 symmetric positive-definite hydrodynamic structure factor matrix 
S}{(k,uj) = (iJU J [5]. The hydrostatic structure factor matrix Sn{k) is obtained by integrating 
Sn(k,uj) over all frequencies, 



S H (k) 



poc T 2 k B T 
Po^bTq 



(13) 



p^c-^bT* 

We use Sn{k) for an ideal gas (i.e., for po = 1, br = 1) to non-dimensionalize Sji(k,uj), for 
example, we express the spatio-temporal cross-correlation between density and velocity through 
the dimensionless hydrodynamic structure factor 



i 



S Pt v(k,uj) = (poc 2 k B To) 2 (p 1 k B T ) 2 (p(k,u)v*(k,u)) . 
For the non-ideal SHSD fluid the density fluctuations have a spectrum 

s-p(fc) = {poc^kBny 1 (p(k)p*(k)) = c^ 2 , 

which only captures the small k behavior of the full (particle) structure factor S(k) (see Fig. [I]), as 
expected of a continuum theory that does not account for the structure of the fluid. Typically only 
the density-density dynamic structure factor is considered because it is accessible experimentally 
via light scattering measurements and thus most familiar. However, in order to fully access the 
validity of the full LLNS system one should examine the dynamic correlations among all pairs of 
variables. The off-diagonal elements of the static structure factor matrix Sff(k) vanish because the 
primitive hydrodynamic variables are instantaneously uncorrelated, however, they have non-trivial 
dynamic correlations visible in the off-diagonal elements of the dynamic structure factor matrix 
S H (k,Lo). 

In Figs. [4] and [5] we compare theoretical and numerical results the hydrodynamic structure 
factors for the SHSD fluid with % = 1 at two densities for a small and a medium k value [kD/(2ir) ~ 
0.01 and 0.08]. In this figure we show selected elements of S}{{k,uj) as predicted by the analytical 



solution to Eqs. (12) with parameters obtained by using the HNC approximation to <?2 in the 



Enskog kinetic theory presented in Section II A2[ Therefore, for SHSD the theoretical calculations 
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Figure 4: Normalized density fluctuations c^S p (k,uj) for /c£> « 0.070 for an ideal Maxwell I-DSMC (4> = 0.5, 
X = 0.62) and two non-ideal SHSD (0 = 0.5, x = 1 and <p = 1, \ = 1) fluids of similar kinematic viscosity, 
as obtained from particle simulations (symbols with parameters , fcsTo — 1, m — !)■ The predictions of 
the LLNS equations are also shown for comparison in the same color (solid lines). For the SHSD fluid we 
obtained the transport coefficients from the Enskog theory with the HNC approximation to 172 , while for the 
Maxwell FDSMC fluid we numerically estimated the viscosity and thermal conductivity. 



of Sn(k,uj) do not use any numerical inputs from the particle runs. We also show hydrodynamic 
structure factors obtained from particle simulations in a quasi-one-dimensional setup in which the 
simulation cell was periodic and long along the x axis, and divided into 60 hydrodynamic cells of 
length 5D. Finite-volume averages of the hydrodynamic conserved variables were then calculated 
for each cell every 10 time steps and a Fast Fourier Transform used to obtain hydrodynamic 
structure factors for several wavenumbers. Figure [4] shows very good agreement between theory 
and numerics, and clearly shows the shifting of the two symmetric Brillouin peaks at u ~ c s k 
toward higher frequencies as the compressibility of the SHSD fluid is reduced and the speed of 
sound increased. Figure [5] shows that the positions and widths of the side Brillouin peaks and the 
width of the central Rayleigh are well-predicted for all elements of Sjj(k,u) for a wide range of k 
values, demonstrating that the SHSD fluid shows the expected fluctuating hydrodynamic behavior. 
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Figure 5: Selected diagonal (left panel) and off-diagonal elements (right panel) of the non-dimensionalizcd 
hydrodynamic structure factor matrix Sn{k,oj) for a large wavenumber kD ss 0.50 for an SHSD fluid at 
(j) = 1, x = 1 (symbols), compared to the predictions from the LLNS equations (lines of same color). The 
remaining parameters are as in Fig. [4] 

C. Brownian Walker VACF 

As an illustration of the correct hydrodynamic behavior of the SHSD fluid and the significance 
of compressibility, we study the velocity autocorrelation function (VACF) C(t) = (v x (0)v x (t)) for 
a single neutrally-buoyant hard sphere Brownian bead of mass M and radius R suspended in an 
SHSD fluid of mass density p. This problem is relevant to the modeling of polymer chains or 
(nano) colloids in solution, and led to the discovery of a long power-law tail in C(t) |38j which has 
since become a standard test for hydrodynamic behavior of solvents \27\ 1391 140] . Here the fluid 
particles interact via stochastic collisions, exactly as in I-DSMC. The interaction between fluid 
particles and the bead is treated as if the SHSD particles are hard spheres of diameter D s , chosen 
to be somewhat smaller than their interaction diameter with other fluid particles (specifically, we 
use D s = D/4) for computational efficiency reasons, using an event-driven algorithm [H]. Upon 
collision with the bead the relative velocity of the fluid particle is reversed in order to provide a 
no-slip condition at the surface of the suspended sphere [iQl HI] (slip boundaries give qualitatively 
identical results). For comparison, an ideal I-DSMC fluid of comparable viscosity is also simulated. 

Theoretically, C(t) has been calculated from the linearized (compressible) fluctuating Navier- 
Stokes (NS) equations [ID] . The results are analytically complex even in the Laplace domain, how- 
ever, at short times an inviscid compressible approximation applies. At large times the compress- 
ibility does not play a role and the incompressible NS equations can be used to predict the long-time 
tail. At short times, t < t c = 2R/c s , the major effect of compressibility is that sound waves gener- 
ated by the motion of the suspended particle carry away a fraction of the momentum, so that the 
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VACF quickly decays from its initial value C(0) = ksT/M to C(t c ) ~ ksl '/M e f f, where M e jj = 
M -\-2itR? p/?>. At long times, t > t v i sc = 4pR H /3r], the VACF decays as in an incompressible fluid, 
with an asymptotic power-law tail (kBT/M)(8\ /r 3ir)~ 1 (t/t v i sc )~ 3 / 2 , in disagreement with predic- 
tions based on the Langevin equation (Brownian dynamics), C(t) = (ksT/M) exp (—6TrRHrjt/M). 
We have estimated the effective (hydrodynamic) colloid radius Rh from numerical measurements 
of the Stokes friction force F = —QitRhtiv and found it to be somewhat larger than the hard-core 
collision radius R + D s /2, but for the calculations below we use Rh = R + D s /2. 




Figure 6: The velocity autocorrelation function for a neutrally buoyant hard sphere suspended in a non-ideal 
SHSD (x = 1) fluid at two densities (symbols), <f> = 0.5 and = 1.0, as well as an ideal Maxwell I-DSMC 
fluid (ip — 0.5, x — 0.62, symbols), at short and long times (inset). For the more compressible (less viscous) 
fluids the long time tails are statistically measurable only up to t/t V i SC ~ 5. The theoretical predictions 
based on the inviscid, for short times, or incompressible, for long times, Navier-Stokes equations are also 
shown (lines). 

In Fig. [6] numerical results for the VACF in a Maxwell I-DSMC fluid and an SHSD fluid at two 
different densities are compared to the theoretical predictions. The diameter of the nano-colloidal 
particle is only 2.5D (i.e., Rh = I.375D), although we have performed simulations using larger 
spheres as well with very similar (but less accurate) results. Since periodic boundary conditions 
were used we only show the tail up to about the time at which sound waves generated by its periodic 
images reach the particle, ti, = L/c s , where the simulation box was L = 25D. In dimensionless 
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units, the viscosity r/ = fjD 2 \JmkBT was measured to be r\ ~ 0.75 for both the Maxwell I-DSMC 
fluids and the SHSD fluid at cfi = 0.5, and fj ~ 1.9 for SHSD at (f> = 1. The results in Fig. [6] are 
averages over 10 runs, each of length T/t visc rs 2 • 10 5 for I-DSMC, T/t visc « 1 • 10 5 for SHSD at 
= 0.5, and T/t v i sc ~ 4.5- 10 4 for SHSD at = 1.0, where in atomistic time units to = D^Jm/ksT 
the viscous time scale is t v i sc /tQ 60/(37177). 

It is seen, as predicted, that the compressibility or the sound speed c s , determines the early 
decay of the VACF. The exponent of the power-law decay at large times is also in agreement with 
the hydrodynamic predictions. The coefficient of the VACF tail agrees reasonably well with the 
hydrodynamic prediction for the less dense fluids, however, there is a significant deviation of the 
coefficient for the densest fluids, perhaps due to ordering of the fluid around the suspended sphere, 
not accounted for in continuum theory. In order to study this discrepancy in further detail one 
would need to perform simulations with a much larger bead. This is prohibitively expensive with 
the serial event-driven algorithm used here jH] and requires either parallelizing the code or using 
a hybrid particle-continuum method [29 , which we leave for future work. 

IV. CONCLUSIONS 

We have successfully generalized the traditional DSMC algorithm for simulating rare gas flows to 
flows of dense non-ideal fluids. Constructing such a thermodynamically-consistent Stochastic Hard 
Sphere Dynamics (SHSD) algorithm required first eliminating the grid artifacts from traditional 
DSMC. These artifacts are small in traditional DSMC simulations of rarefied gases because the col- 
lisional cell size is kept significantly smaller than the mean free path |42 j , but become pronounced 
when dense flows are simulated because the collisional-stress tensor is not isotropic. Our Isotropic 
DSMC (I-DSMC) method is a grid free DSMC variant with pairwise spherically-symmetric stochas- 
tic interactions between the particles, just as classical fluids simulated by molecular dynamics (MD) 
use a pairwise spherically-symmetric deterministic interaction potential. The I-DSMC method can 
therefore be viewed as a transition from the DSMC method, suitable for rarefied flows, to the MD 
method, suitable for simulating dense liquids (and solids). 

It has long been apparent that manipulating the stochastic collision rules in DSMC can lead 
to a wide range of fluid models, including non-ideal ones |20| |4"3] . It has also been realized that 
DSMC, as a kinetic Monte Carlo method, is not limited to solving the Boltzmann equation [31] 
but can be generalized to Enskog-like kinetic equations \22\ [2"3"] . However, what has been so far 
elusive is to construct a DSMC collision model that is thermodynamically consistent, meaning that 
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the resulting fluid structure and the equation of state are consistent with each other as required 
by statistical mechanics. We overcame this hurdle here by constructing stochastic collision kernels 
in I-DSMC to be as close as possible to those of the classical hard-sphere deterministic system. 
Thus, in the SHSD algorithm randomly chosen pairs of approaching and overlapping particles 
undergo collisions as if they were hard spheres of variable diameter. This is similar to the modified 
collision rules used to construct a consistent non-ideal Multi Particle Collision Dynamics fluid in 
Refs. [I811I9]. 

We demonstrated the consistent thermodynamic behavior of the SHSD system by observing that 
it has identical structure and thermodynamic properties to a Hamiltonian system of penetrable 
spheres interacting with a linear core potential, even up to solid densities. We found that at 
fluid densities the pair correlation function 52 (f) of the linear core system is well-described by the 
approximate HNC closure, enabling us to obtain moments of 52 (?")■ These moments were then used 
as inputs in a modified Chapman-Enskog calculation to obtain excellent estimates of the equation 
of state and transport coefficients of the SHSD fluid over a wide range of densities. We do not yet 
have a complete theoretical understanding of our surprising finding that the SHSD system behaves 
thermodynamically identically to the linear core system. An important open question remains 
whether by choosing a different collision kernel one can obtain stochastic fluids corresponding to 
Hamiltonian systems of penetrable spheres interacting with effective pair potentials U e ff{r) other 
than the linear core potential. 

The SHSD algorithm is similar in nature to DPD and has a similar computational complexity. 
The essential difference is that DPD has a continuous-time formulation (a system of stochastic 
ODEs), where as the SHSD dynamics is discontinuous in time (Master Equation). This is similar 
to the difference between MD for continuous potentials and discontinuous potentials. Just as 
DSMC is a stochastic alternative to hard-sphere MD for low-density gases, SHSD is a stochastic 
modification of hard-sphere MD for dense gases. On the other hand, DPD is a modification of MD 
for smooth potentials to allow for larger time-steps and a conservative thermostat. 

A limitation of SHSD is that for reasonable values of the collision frequency (x ~ 1) and density 
(cf) ~ 1) the fluid is still relatively compressible compared to a dense liquid, S(k = 0) = c^ 2 > 0.1. 
Indicative of this is that the diffusion coefficient is large relative to the viscosity as it is in typical 
DPD simulations, so that the Schmidt number S c = ^(pC) 1 is l ess than 10 instead of being on the 
order of 100-1000. Achieving higher ct or S c requires high collision rates (for example, % ~ 10 4 
is used in Ref. [32 ) and appropriately smaller time steps to ensure that there is at most one 
collision per particle per time step, and this requires a similar computational effort as in hard- 
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sphere molecular dynamics at a comparable density. At low and moderate gas densities the SHSD 
algorithm is not as efficient as DSMC at a comparable collision rate. However, for a wide range 
of compressibilities, SHSD is several times faster than the alternative deterministic Event-Driven 
MD (EDMD) for hard spheres |24[ I44j . Furthermore, SHSD has several important advantages over 
EDMD, in addition to its simplicity: 

1. SHSD has several controllable parameters that can be used to change the transport coeffi- 
cients and compressibility, notably the usual density (f> but also the cross-section factor x 
and others 3 , while EDMD only has density. 

2. SHSD is time-driven rather than event-driven thus allowing for easy parallelization. 

3. SHSD can be more easily coupled to continuum hydrodynamic solvers, just like ideal-gas 
DSMC [U] and DPD [461 147] . Strongly-structured particle systems, such as fluids with strong 
interparticle repulsion (e.g., hard spheres), are more difficult to couple to hydrodynamic 
solvers [IB] than ideal fluids, such as MPCD or (I-)DSMC, or weakly-structured fluids, such 
as DPD or SHSD fluids. 

Finally, the stochastic particle model on which SHSD is based is intrisically interesting and the- 
oretical results for models of this type will be helpful for the development of consistent particle 
methods for fluctuating hydrodynamics. 
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